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<N ! ABSTRACT 

s : 

This paper addresses the problem of giving conditions for transcriptional systems to be glob- 
ally entrained to external periodic inputs. By using contraction theory, a powerful tool from 
dynamical systems theory, it is shown that certain systems driven by external periodic signals 
have the property that all solutions converge to a fixed limit cycle. General results are proved, 
and the properties are verified in the specific case of some models of transcriptional systems. 

o\ 

O ; 1 Introduction 

x> : 

Periodic, clock-like rhythms pervade nature and regulate the function of all living organisms. 

For instance, circadian rhythms are regulated by an endogenous biological clock entrained by 
, the light signals from the environment that then acts as a pacemaker, (Gonze et al. 2005). 

Moreover, such an entrainment can be obtained even if daily variations are present, like e.g. 

temperature and light variations. Another important example of entrainment in biological 
, systems is at the molecular level, where the synchronization of several cellular processes is 

regulated by the cell cycle (Tyson et al. 2002). 

An important question in mathematical and computational biology is that of finding condi- 
Q\ | tions ensuring that entrainment occurs. The objective is to identify classes of biological systems 

that can be entrained by an exogenous signal. To solve this problem, modelers often resort to 
£> . simulations in order to show the existence of periodic solutions in the system of interest. Simu- 

lations, however, can never prove that solutions will exist for all parameter values, and they are 
subject to numerical errors. Moreover, robustness of entrained solutions needs to be checked in 
the presence of noise and uncertainties, which cannot be avoided experimentally. 

From a mathematical viewpoint, the problem of formally showing that entrainment takes 
place is known to be very difficult. Indeed, if a stable linear time-invariant model is used to 
represent the system of interest, then entrainment is usually expected, when the system is driven 
by an external periodic input, with the system response being a filtered, shifted version of the 
external driving signal. However, in general, as is often the case in biology, models are nonlinear. 
The response of nonlinear systems to periodic inputs is the subject of much current systems 
biology experimentation; for example, in (Mettetal et al. 2008), the case of a cell signaling 
system driven by a periodic square-wave input is considered. From measurements of a periodic 
output, the authors fit a transfer function to the system, implicitly modeling the system as 
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linear even though (as stated in the Suppemental Materials to (Mettetal et al. 2008)) there are 
saturation effects so the true system is nonlinear. For nonlinear systems, driving the system 
by an external periodic signal does not guarantee the system response to also be a periodic 
solution, as nonlinear systems can exhibit harmonic generation or suppression and complex 
behaviour such as chaos or quasiperiodic solutions (Kuznetsov 2004). This may happen even 
if the system is well-behaved with respect to constant inputs; for example, there are systems 
which converge to a fixed steady state no matter what is the input excitation, so long as this 
input signal is constant, yet respond chaotically to the simplest oscillatory input; we outline 
such an example in an Appendix to this paper, see also (Sontag 2009). Thus, a most interesting 
open problem is that of finding conditions for the entrainment to external inputs of biological 
systems modelled by sets of nonlinear differential equations. 

One approach to analyzing the convergence behavior of nonlinear dynamical systems is 
to use Lyapunov functions. However, in biological applications, the appropriate Lyapunov 
functions are not always easy to find and, moreover, convergence is not guaranteed in general 
in the presence of noise and/or uncertainties. Moreover, such an approach can be hard to 
apply to the case of non-autonomous systems (that is, dynamical systems directly dependent 
on time), as is the case when dealing with periodically forced systems. 

The above limitations can be overcome if the convergence problem is interpreted as a 
property of all trajectories, asking that all solutions converge towards one another (contrac- 
tion). This is the viewpoint of contraction theory, (Lohmiller & Slotine 1998), (Lohmiller & 
Slotine 2000), and more generally incremental stability methods (Angeli 2002). Global results 
are possible, and these are robust to noise, in the sense that, if a system satisfies a contraction 
property then trajectories remain bounded in the phase space (Pham et al. 2009). Contraction 
theory has a long history. Contractions in metric functional spaces can be traced back to the 
work of Banach and Caccioppoli (Granas &: Dugundji 2003) and, in the field of dynamical sys- 
tems, to (Hartman 1961) and even to (Lewis 1949) (see also (Pavlov et al. 2004), (Angeli 2002), 
and e.g. (Lohmiller & Slotine 2005) for a more exhaustive list of related references). Contrac- 
tion theory has been successfully applied to both nonlinear control and observer problems, 
(Lohmiller & Slotine 2000), (Jouffroy &: Slotine 2004) and, more recently, to synchronization 
and consensus problems in complex networks (Slotine et al. 2004), (Wang & Slotine 2005). 
In (Russo and di Bernardo 2009) it was proposed that contraction can be particularly useful 
when dealing with the analysis and characterization of biological networks. In particular, it 
was found that using non Euclidean norms can be particularly effective in this context (Russo 
and di Bernardo 2009), (Russo and di Bernardo 2009b). 

One of the objectives of this paper is to give a self-contained exposition, with all proofs 
included, of results in contraction theory as applied to entrainment of periodic signals, and, 
moreover, to show their applicability to a problem of biological interest, having to do with a 
driven transcriptional system. A surprising fact is that, for these applications, and contrary 
to many engineering aplications, norms other than Euclidean, and associated matrix measures, 
must be considered. 

1.1 Mathematical tools 

We consider in this paper systems of ordinary differential equations, generally time-dependent: 

x = f(t,x) (1) 
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defined for t £ [0, oo) and x S C, where C is a subset of MJ 1 . It will be assumed that f(t,x) is 
differentiable on x, and that f(t,x), as well as the Jacobian of / with respect to x, denoted as 
J(t, x) = x), are both continuous in (i, x). In applications of the theory, it is often the case 
that C will be a closed set, for example given by non-negativity constraints on variables as well 
as linear equalities representing mass-conservation laws. For a non-open set C, differentiability 
in x means that the vector field f(t, •) can be extended as a differentiable function to some 
open set which includes C, and the continuity hypotheses with respect to (t, x) hold on this 
open set. 

We denote by (p(t, s, £) the value of the solution x(t) at time t of the differential equation ((!]) 
with initial value x(s) = £. It is implicit in the notation that (f(t, s,£) € C ("forward invariance" 
of the state set C). This solution is in principle defined only on some interval s < t < s + e, 
but we will assume that <p(t,s,£) is defined for all t > s. Conditions which guarantee such 
a "forward-completeness" property are often satisfied in biological applications, for example 
whenever the set C is closed and bounded, or whenever the vector field / is bounded. (See 
Appendix C in (Sontag 1998) for more discussion, as well as (Angeli & Sontag 1999) for a 
characterization of the forward completeness property.) Under the stated assumptions, the 
function tp is jointly differentiable in all its arguments (this is a standard fact on well-posedness 
of differential equations, see for example Apendix C in (Sontag 1998)). 

We recall (see for instance (Michel et al. 2007)) that, given a vector norm on Euclidean 
space ([•!), with its induced matrix norm \\A\\, the associated matrix measure fi is defined as 
the directional derivative of the matrix norm, that is, 

==jtas(IIJ + M|-i). 

For example, if |»| is the standard Euclidean 2-norm, then fi(A) is the maximum eigenvalue of 
the symmetric part of A. As we shall see, however, different norms will be useful for our appli- 
cations. Matrix measures are also known as "logarithmic norms" , a concept independently in- 
troduced by Germund Dahlquist and Sergei Lozinskii in 1959, (Dahlquist 1959, Lozinskii 1959). 
The limit is known to exist, and the convergence is monotonic, see (Strom 1975, Dahlquist 1959). 

We will say that system ([TJ is infinitesimally contracting on a convex set C C R n if there 
exists some norm in C, with associated matrix measure [i such that, for some constant c £ 
M-{0}, 

fi(J(x,t))<-c 2 , VxeC, W>0. (2) 

Let us discuss very informally (rigorous proofs are given later) the motivation for this con- 
cept. Since by assumption / (t, x) is continuously differentiable, the following exact differential 
relation can be obtained from ([1]): 

5x = J (t, x) 5x, (3) 

where, as before, J = J (t, x) denotes the Jacobian of the vector field /, as a function of x € C 
and t € M + . (The object 5x can be thought of as a "virtual displacement" in the language of 
mechanics, as in (Arnold 1978), which views such displacements as linear tangent differential 
forms differentiable with respect to time.) Consider now two neighboring trajectories of ([T]), 
evolving in C, and the virtual displacements between them. Note that ([3]) can be seen as a 
linear time- varying dynamical system of the form: 

5x = J (t) 6x. 
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Hence, an upper bound for the magnitude of its solutions can be obtained by means of the 
Coppel inequality (Vidyasagar 1993), yielding: 

\6x\<\6x \efirt J ®W, (4) 

where \x (J) is the matrix measure of the system Jacobian induced by the norm being considered 
on the states and \5x (0)| = \Sxo\. Using ^ and ([2]), we have that 

3 /3>0: \6x(t)\ < pe~ c2t . 

Thus, trajectories starting from infinitesimally close initial conditions converge exponentially 
towards each other. In what follows we will refer to c 2 as contraction (or convergence) rate. 

The key theoretical result about contracting systems links infinitesimal and global contrac- 
tivity, and is stated below. This result can be traced, under different technical assumptions, to 
e.g. (Lohmiller & Slotine 1998), (Pavlov et al. 2004), (Lewis 1949), (Hartman 1961). 

Theorem 1. Suppose that C is a convex subset ofM. n and that f(t, x) is infinitesimally contract- 
ing with contraction rate c 2 . Then, for every two solutions x(t) = ip(t, 0, £) and z(t) = c/?(i, 0, () 
of (G]), it holds that: 

\x(t)-z{t)\ < e" c2 *|e-CI, Vi>0. (5) 

In other words, infinitesimal contr activity implies global contr activity. In the Appendix, we 
provide a self-contained proof of Theorem [TJ In fact, the result is shown there in a generalized 
form, in which convexity is replaced by a weaker constraint on the geometry of the space. 

In actual applications, often one is given a system which depends implicitly on the time, t, by 
means of a continuous function u (t), i.e. systems dynamics are represented by x = f (x, u (t)). 
In this case, u (t) : M + — > U (where U is some subset of K), represents an external input. It 
is important to observe that the contractivity property does not require any prior information 
about this external input. In fact, since u (t) does not depend on the system state variables, 
when checking the property, it may be viewed as a constant parameter, u € U. Thus, if 
contractivity of / (x, u) holds uniformly Vit £ U, then it will also hold for / (x, u (t)). 

Given a number T > 0, we will say that system ([TJ is T -periodic if it holds that 

f(t + T,x) = f(t,x) Vt>0, xeC. 

Notice that the system x = f(x,u(t)) is T-periodic, if the external input, u(t), is itself a 
periodic function of period T. 

The following is the basic theoretical result about periodic orbits that will be used in the 
paper. It may be found, under various different technical variants, in the references given above. 

Theorem 2. Suppose that: 

• C is a closed convex subset ofM. n ; 

• f is infinitesimally contracting with contraction rate c 2 ; 

• f is T-periodic. 

Then, there is a unique periodic solution a(t) : [0, oo) — > C of ([!]) of period T and, for every 
solution x(t), it holds that \x (t) — a (t)\ — > as t — > oo. 

In the Appendix of this paper, we provide a self-contained proof of Theorem in a gener- 
alized form which does not require convexity. 
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1.2 A simple example 



As a first example to illustrate the application of the concepts introduced so far, we choose a 
simple bimolecular reaction, in which a molecule of A and one of B can reversibly combine to 
produce a molecule of C. 

This system can be modeled by the following set of differential equations: 

A= -kyAB + k-iC, 

B = -k x AB + k-iC, (6) 
C = k\AB - k-iC, 

where we are using A = A(t) to denote the concentration of A and so forth. The system evolves 
in the positive orthant of R 3 . Solutions satisfy (stoichiometry) constraints: 

A(t) + C{t) = a 

B(t) + C(t) = P {) 

for some constants a and (3. 

We will assume that one or both of the "kinetic constants" k{ are time-varying, with period 
T. Such a situation arises when the k^s depend on concentrations of additional enzymes, 
which are available in large amounts compared to the concentrations of A, B, C, but whose 
concentrations are periodically varying. The only assumption will be that ki(t) > k® > and 
k-i(t) > k\ > for all t. 

Because of the conservation laws ([7]), we may restrict our study to the equation for C. Once 
that all solutions of this equation are shown to globally converge to a periodic orbit, the same 
will follow for A(t) =a- C(t) and B(t) = (3 - C{t). We have that: 

C = h(a-C)(P -C) -fc_iC. (8) 

Because A(t) > and B(t) > 0, this system is studied on the subset of R defined by < C < 
min{a,/5}. The equation can be rewritten as: 

C = k 1 (ap-aC-pC + C 2 ) -k^C. (9) 

Differentiation with respect to C of the right-hand side in the above system yields this (lxl) 
Jacobian: 

J := An (- (a + P) + 2C - k-i) . (10) 
Since we know that —a + C < and — j3 + C < 0, it follows that 

J < -Jfeifc-i < -Jfejife" ! := -c 2 

for c = ■J k^k^i- Using any norm (this example is in dimension one) we have that /i( J) < — c 2 . 
So ([6]) is contracting and, by means of Theorem [21 solutions will globally converge to a unique 
solution of period T (notice that such a solution depends on system parameters). 

Figure [T] shows the behavior of the dynamical system Q, using two different values of k—\. 
Notice that the asymptotic behavior of the system depends on the particular choice of the 
biochemical parameters being used. Furthermore, it is worth noticing here that the higher the 
value of k-i, the faster will be the convergence to the attractor. 
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Figure 1: Entrainment of © to the periodic input u(t) = 1.5 + sin(lOi) for k-i = 10 (blue), 
k-i = 1 (green), k-i = 0.1 (red). Notice that an increase of k-i, causes an increase of the 
contraction rate, hence trajectories converge faster to the system unique periodic attractor. 
The other system parameters are set to: a = (3 = 1, k-2 = 0.1. 

2 Results 

2.1 Mathematical model and problem statement 

We study a general externally-driven transcriptional module. We assume that the rate of 
production of a transcription factor X is proportional to the value of a time dependent input 
function u(t), and X is subject to degradation and/or dilution at a linear rate. (Later, we 
generalize the model to also allow nonlinear degradation as well.) The signal u(t) might be an 
external input, or it might represent the concentration of an enzyme or of a second messenger 
that activates X. In turn, X drives a downstream transcriptional module by binding to a 
promoter (or substrate), denoted by E, whose free concentration is denoted as e = e(t). The 
binding reaction of X with E is reversible and given by: 



where Y is the complex protein-promoter, and the binding and dissociation rates are k\ and £?2 
respectively. As the promoter is not subject to decay, its total concentration, Et, is conserved, 
so that the following conservation relation holds: 



We wish to study the behavior of solutions of the system that couples X and E, and specifically 
to show that, when the input u(t) is periodic with period T, this coupled system has the property 
that all solutions converge to some globally attracting limit cycle whose period is also T. 

Such transcriptional modules are ubiquitous in biology, natural as well as synthetic, and 
their behavior was recently studied in (Del Vecchio et al. 2008) in the context of "retroactivity" 
(impedance or load) effects. If we think of u(t) as the concentration of a protein Z that 
is a transcription factor for X, and we ignore fast mRNA dynamics, such a system can be 
schematically represented as in Figure[2l which is adapted from (Del Vecchio et al. 2008). Notice 



X + E ^ Y, 



E + Y = E T . 



(11) 
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Figure 2: A schematic diagram of the two transcriptional modules modeled in (|12p 

that u(t) here does not need to be the concentration of a transcriptional activator of X for our 
results to hold. The results will be valid for any mathematical model for the concentrations, x, 
of X and y, of Y (the concentration of E is conserved) of the form: 

x = u (t) - 5x + kiy - k 2 (E T - y) x . 
y = -hy + k 2 {E T -y)x. 

Our main objective in this paper is, thus, to show that, when u is a periodic input, all 
solutions of system (j!2jl converge to a (unique) limit cycle (Figure [3]) . The key tool in this 
analysis is to show that, when no input is present, the system is infinitesimally, and hence 
globally, contracting. 

Thus, the main step will be to establish the following technical result, see Section [2.21 
Theorem 3. The system 

x = — 8x + k\y — k 2 (Et — y) x 
y = -kiy + k 2 (E T -y)x 

where 

(x(t),y(t))eC = [0,oo)x[0,E T ] (13) 
for all t > 0, and Ex, k±, k 2 , and 5 are arbitrary positive constants, is contracting. 

By means of Theorem [2j we then have the following immediate Corollary: 

Theorem 4. For any given nonnegative periodic input u of period T, all solutions of system 
Ql| ) converge exponentially to a periodic solution of period T. 



In the following sections, we introduce a matrix measure that will help establish contrac- 
tivity, and we prove Theorem [3l We will also discuss several extensions of this result, allowing 
the consideration of multiple driven subsystems as well as more general nonlinear systems with 
a similar structure. 



2.2 Proof of Theorem H 



We will use Theorem [2j The Jacobian matrix to be studied is: 

-5 - k 2 (Et -y) ki + k 2 x 
k 2 (E T -y) -h- k 2 x 



J := 



(14) 
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Figure 3: Entrainment of the transcriptional module (|12p output (green), Y, to the periodic 
input (blue): u(t) = 1.5 + sin(O.lt) (left) and to a repeating {0,1} sequence (right). System 
parameters are set to: 5 = 3, k%=l, k<i = 0.1. 
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As matrix measure, we will use the measure /x^i induced by the vector norm [Psc^, where P 
is a suitable nonsingular matrix. More specifically, we will pick P diagonal: 



Pi 
p 2 



(15) 



where p\ and p 2 are two positive numbers to be appropriately chosen depending on the param- 
eters defining the system. 

It follows from general facts about matrix norms that 



fx P ,x (J) = n\ {PJP 



-V 



(16) 



where [i\ is the measure associated to the 1*^ norm and is explicitly given by the following 
formula: 



Mi (J) = max Jjj + ^2 \Jij\ 

&3 



(17) 



Observe that, if the entries of J are negative, then asking that /Ui(J) < amounts to a column 
diagonal dominance condition. (The above formula is for real matrices. If complex matrices 
would be considered, then the term Jjj should be replaced by its real part $l{Jjj}.) 

Thus, the first step in computing \xp \ (J) is to calculate PJP^ 1 : 



% [hi (E T - y)] -h - k 2 x 



pi 



(18) 



Using (fTTj) . we obtain: 



fJ>p,l (J) = max 1—5 - k 2 {E T — y) + 



— k 2 (Et - y) 
Pi 



; —ki - k 2 x + 



— (ki + k 2 x) 

P2 



(19) 



Note that we are not interested in calculating the exact value for the above measure, but just 
in ensuring that it is negative. To guarantee that \xp^\ (J) < 0, the following two conditions 
must hold: 



5 - k 2 {Et -y) + 



—k 2 (E T - y) 
Pi 



< -q; 



ki - k 2 x + 



— (ki + k 2 x) 

P2 



< ~<?2 



(20) 



(21) 



Thus, the problem becomes that of checking if there exists an appropriate range of values for 
pi, p 2 that satisfy (f20|) and (f2TI) simultaneously. 

The left hand side of (|21|) can be written as: 



P±-l) (ki + k 2 x), 

P2 J 

which is negative if and only if p\ < p 2 . In particular, in this case we have: 



(22) 



(kx + k 2 x) <l--l)h:-- 

P2 ) \P2 
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The idea is now to ensure negativity of (|20p by using appropriate values for p\ and p 2 which 
fulfill the above constraint. Recall that the term Et — y > because of the choice of the state 
space (this quantity represents a concentration). Thus, the left hand side of ([2D]) becomes 

-6+(^-l\k 2 (E T -y) (23) 

The next step is to choose appropriately p 2 and p\ (without violating the constraint p 2 > pi). 
Imposing p 2 /p\ = 1 + £, £ > 0, (|2"5j) becomes 

- 5 + ek 2 (E T - y) . (24) 

Then, we have to choose an appropriate value for e in order to make the above quantity 
uniformly negative. In particular, (|24p is uniformly negative if and only if 

5 5 

£< k 2 (E T -y) -T^Et (25) 

We can now choose 

5 



k 2 E T 

with < £ < . In this case, (|24|l becomes 

-5 + ek 2 (E T - y) < -£k 2 E T := -c|. 

Thus, choosing p\ = 1 andj>2 = 1+e = 1 + k ^ Er with < £ < k ^ ET ' we nave (J) < — ° 2 - 
Furthermore, the contraction rate c 2 , is given by: 



mm 



Rc 2 .} 



Notice that c 2 depends on both system parameters and on the elements pi, p 2 , i.e. it depends on 
the particular metric chosen to prove contraction. This completes the proof of the Theorem. □ 

2.3 Generalizations 

In this Section, we discuss various generalizations that use the same proof technique. 
2.3.1 Assuming X activation by enzyme kinetics 

The previous model assumed that X was created in proportion to the amount of external signal 
u(t). While this may be a natural assumption if u(t) is a transcription factor that controls the 
expression of X, a different model applies if, instead, the "active" form X is obtained from 
an "inactive" form Xq, for example through a phosphorylation reaction which is catalyzed by 
a kinase whose abundance is represented by u(t). Suppose that X can also be constitutively 
deactivated. Thus, the complete system of reactions consists of 

X + E ^ Y, 

together with 

X n ^X 
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where the forward reaction depends on u. Since the concentrations of X$ + X + Y must remain 
constant, let us say at a value ^tot> we eliminate Xq and have: 

x = u(t)(X to t - x - y) - 5x + fay - k 2 (E T - y)x, , , 

V = -fay + k 2 (E T -y)x. 

We will prove that if u (t) is periodic and positive, i.e. u (t) > u$ > 0, then a globally 
attracting limit cycle exists. Namely, it will be shown, after having performed a linear coordinate 
transformation, that there exists a negative matrix measure for the system of interest. 

Consider, indeed, the following change of the state variables: 

x t = x + y. (27) 

The systems dynamics, then become: 

x t = u (t) (X t ot - x t ) - 5x t + 5y , , 

y = ~hy + k 2 (E T - y) (x t - y) 

As matrix measure, we will now use the measure /Uqo induced by the vector norm 1*1^. (Notice 
that this time, the matrix P is the identity matrix). 

Given a real matrix J, the matrix measure /ioo (J) is explicitly given by the following formula 
(see e.g. (Michel et al. 2007)): 

(J) = max \Jn + ^2\Jij\\ ■ (29) 
V & ) 

(Observe that this is a row-dominance condition, in contrast to the dual column-dominance 
condition used for 

Differentiation of (|28p yields the Jacobian matrix: 

-u(t)-S 5 
k 2 (Et - y) -fa + k 2 (-E T -x t + 2y) 



J :-- 



Thus, it immediately follow from (|29l) that fi^ (J) is negative if and only if: 

-u(t)-8+\S\<-4; (30) 

- fa + k 2 {-E T -x t + 2y) + \k 2 (Et - y)\ < -c|. (31) 

The first inequality is clearly satisfied since by hypotheses both system parameters and the 
periodic input u (t) are positive. In particular, we have: 

— u (t) — 5 + \5\ < —uq := — cf ; 

By using ([271) (recall that Et — y > 0), the right hand side of the second inequality can be 
written as: 

— fa + k 2 (-E T - x t + 2y) + k 2 (E T — y) = - fa - fax. 

Since all system parameters are positive and x > 0, the above quantity is negative and upper 
bounded by —fa := — c|. 
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Thus, W6 have that /ioo 

(J) < — c 2 , where: 

2 r 2 21 

c = mm ( Ci , c 2 j . 

The contraction property for the system is then proved. By means of Theorem [2j we can then 
conclude that the system can be entrained by any periodic input. 

Simulation results are presented in Figure HI where the presence of a stable limit cycle 
having the same period as u (t) is shown. 




50 100 150 

time (minutes) 




01 23456789 10 

time (minutes) 



Figure 4: Left panel: entrainment of the transcriptional module (|26p output (green), Y, to the 
periodic input (blue): u(t) = 1.5 + sin(O.lt). Right panel: zoom on t S [0,10] min. showing 
trajectories starting from different initial conditions converging towards the attracting limit 
cycle. System parameters are set to: k% = 0.5, £?2 = 5, X to t = 1, Ex = 1, 6 = 20. 



2.3.2 Multiple driven systems 

We may also treat the case in which the species X regulates multiple downstream transcrip- 
tional modules which act independently from each other, as shown in Figure The biochemical 
parameters defining the different downstream modules may be different from each other, repre- 
senting a situation in which the transcription factor X regulates different species. After proving 
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Figure 5: A schematic diagram of the transcriptional modules given in (|12|) 

a general result on oscillations, and assuming that parameters satisfy the retroactivity estimates 
discussed in (Del Vecchio et al. 2008), one may in this fashion design a single input-multi output 
module in which e.g. the outputs are periodic functions with different mean values, settling 
times, and so forth. 

We denote by E%, . . . , E n the various promoters, and use yi, . . . ,y n to denote the concen- 
trations of the respective promoters complexed with X. The resulting mathematical model 
becomes: 

x = u(t) - Sx + K n yx - K 2i(-Er,i - yi)x + 

+ K 12 y 2 - K 2 2{E T ,2 ~ V2)x + ■■■ 

+ Ki n y n - K 2n (ET, n - y n )x 
yi = -K uyi + K 2 i(E Tjl - yi )x ( 32 ) 



Vn 



-K ln y n + K 2n (E T ,n ~ Vn)X . 



We consider the corresponding system with no input first, assuming that the states satisfy 
x(t) > and < yi(t) < E T)i for all t,i. 

Our generalization can be stated as follows: 

Theorem 5. System with no input (i.e. u(t) = 0) is contracting. Hence, ifu(t) is a non- 
zero periodic input, its solutions exponentially converge towards a periodic orbit of the same 
period as u(t). 

Proof. We only outline the proof, since it is similar to the proof of Theorem We employ the 
following matrix measure: 

fi PA (J) = /i! (PJP- 1 ) , (33) 
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where 



P :-- 



pi ... 
p 2 ... 



... p n +i 
and the scalars pi have to be chosen appropriately (pi > 0, Vi = 1, . . 



In this case, 



J : = 



and 



K 2 i(E Ttl - yi) 
K 2 2(Et,2 — IJ2) 

,n Vn) 



Ku + K21X K12 + K22X 
-K u - K21X 

-K 12 - K 22 x 



,n+l). 

Kin + K 2n x 



-Ki n ~ K 2n x 



(34) 



(35) 



PJP- 1 



^2l(^T,l-2/l) 
^K 22 (Et,2 ~ 2/2) 



pi 



p2 (K n +K 2 ix) f 3 (K l2 +K22x) 



-Ku - K21X 






-iTi2 - -^222; 



Pn + l 
Pi 



K 2n (Ex 



Hence, the n + 1 inequalities to be satisfied are: 



n 1 " 

6 - V if2i(^T,i - ^) + — |#2i(#T. 
i=l ■ Pi i=l 



and 



- Ku - K 2 iX + 



Pi 



P1+1 



(K u + K 2i )x 



< -c 



■i+l> 



yOl < -ci 



i = 1, 2, . . . , n. 






-Kl„ - i^2n^ 
(36) 



(37) 



(38) 



Clearly, the set of inequalities above admits a solution. Indeed, the left hand side of ([38 
can be recast as 

' Pi 



P1+1 



1 {Ku + K 2l x) 



i = 1,2 



n 



which is negative definite if and only if pi/pi + i < 1 for all i = 1, . . . , n. Specifically, in this case 
we have 



■^--l) {K li + K 2i x)< [^-\ \ K 



i = 1,2, 



,n 



,Pi+l J \Pi+l 
Also, from ([37]) . as Et % — Vi > for all i, we have that ([37]) can be rewritten as: 



-6 - T K^Er.i - Vi ) + V ^(E T4 - Vi ) < -c\. 
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Since Pi/pi+i < 1, we can impose Pi+i/pi = l+e^j+i (with £1,1+1 > 0) and the above inequality 
becomes 

n 

-5 +} j ei > i +1 K 2 i(E Tt i - yi) < -c\. 

i=l 

Clearly, such inequality is satisfied if we choose £i %+i sufficiently small; namely: 

5 



(n — 1) k 2 Ex,i 

Following a similar derivation to that of Section 12.21 we can choose 

S 

£i+1 " (n - 1) k 2 E T , 4m ' 
with < < ( n _ 1 - ) fc 2 Et ■ In this case, we have: 



' n — 1 



j=l 



Thus, /i(J) < — c 2 , where 

c 2 = min {q} , i = 1, . . . , n + 1. 

i 

The second part of the theorem is then proven by applying Theorem [2j □ 

In Figure[6]the behavior of two-driven downstream transcriptional modules is shown. Notice 
that both the downstream modules are entrained by the periodic input u(t), but their steady 
state behavior is different. 

Notice that, by the same arguments used above, it can be proven that 

x = u{t) (X TO t -x- Ya=i Vi) ~Sx + Knyi - K 2 i(E T ,i - yi)x + 
+ K 12 y 2 - K 22 (E Tt2 - y 2 )x + • • • 
+ Ki n y n - K 2n (E T 

,n yn)X 

yi = -Kim + K 21 {E Tj i - yi)x ' ( 39 ) 

y n = -K ln y n + K 2n (E Tjn - y n )x . 
is contracting. 

2.3.3 Transcriptional cascades 

A cascade of (infinitesimally) contracting systems is also (infinitesimaHy) contracting (see Ap- 
pendix [D] for the proof). This implies that any transcriptional cascade, will also give rise to 
a contracting system, and, in particular, will entrain to periodic inputs. By a transcriptional 
cascade we mean a system as shown in Figure [3 In this figure, we interpret the intermediate 
variables Yi as transcription factors, making the simplifying assumption that TF concentration 
is proportional to active promoter for the corresponding gene. (More complex models, incor- 
porating transcription, translation, and post-translational modifications could themselves, in 
turn, be modeled as cascades of contracting systems.) 
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50 100 150 

time (minutes) 



0.18 




0.02 



time (minutes) 

Figure 6: Outputs Y\ and Yi of two transcriptional modules driven by the external periodic 
input u(t) = 1.5 + sin(i). The parameters are set to: 5 = 0.01, kn = 10, k^x = 10, Et,i = 1 for 
module 1 and k±2 = 0.1, /C22 =0.1, Et,2 = 1 for module 2. 




Figure 7: Transcriptional cascade discussed in the text. Each red box contains the transcrip- 
tional module described by (|12p 
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2.3.4 More abstract systems 



We can extend our results even further, to a larger class of nonlinear systems, as long as the 
same general structure is present. This can be useful for example to design new synthetic 
transcription modules or to analyze the entrainment properties of general biological systems. 
We start with a discussion of a two dimensional system of the form: 



y 



u(t) -a(x) + f(x,y), 
-f(x,y) ■ 



(40) 



In molecular biology, a(x) would typically represent a nonlinear degradation, for instance in 
Michaelis-Menten form, while the function / represents the interaction between x and y. The 
aim of this Section is to find conditions on the degradation and interaction terms that allow 
one to show contractivity of the unforced (no input u) system, and hence existence of globally 
attracting limit cycles. 

We assume that the state space C is compact (closed and bounded) as well as convex. 

Theorem 6. System ( f^0| ), without inputs u, evolving on a convex compact subset of phase 
space is contracting, provided that the following conditions are all satisfied, for each x,y G C: 



• t > 0; 

. 9f 



does not change sign; 



da o<9/ 
dx > z dx ■ 



Notice that the last condition is automatically satisfied if ^ < 0, because > 0. 



(I.r 



Proof. As before, we prove contraction by constructing an appropriate negative measure for 
the Jacobian of the vector field. In this case, the Jacobian matrix is: 



J 



.da + 91 9£ 

dx dx dv 

_£/ _££ 

dx dy 



Once again, as matrix measure we will use: 

A*P,l (J) = Ml {PJP' 1 ) , 

with 



Pi 
p 2 



(41) 

(42) 
(43) 



and pi,p2 > appropriately chosen. 
Using (|4*2|) we have 



, da df 
MP,iO/)=maxi-- + - + 



P2 df_ 
pi dx 



df 



dy 



+ 



Pidf_ 
p 2 dy 



(44) 
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Following the same steps as the proof of Theorem [31 we have to show that: 



- a -i + 

dy 



Pidf_ 
P2 dy 



< -4, 



da df 
h — + 

dx dx 



P2,df_ 
Pi dx 



< -4 



(45) 



(46) 



Clearly, if df /dy > for every x,y G C and p\ < P2, the first inequality is satisfied, with 



P2 J dx' 



To prove the theorem we need to show that there exists p\ < p2 and c 2 . satisfying (|46p . For 
such inequality, since df /dx does not change sign in C by hypothesis, we have two possibilities: 

1. g <0, Vx,y£C; 

2. g > , Vx,yGC. 

In the first case, the right hand side of (|46p becomes 



(47) 



dx dx p\ dx 



Choosing P2/P1 = 1 + e, with e > 0, we have: 



da df_ _P2_df_ _ _da + £ df_ 
dx dx pi dx dx dx 



Specifically, if we now pick 



where A = max |^ and B = min 
definite, i.e. 



e > 



A 
~B 



, we have that the above quantity is uniformly negative 



2 da df 2 

dC2 < i: ~di + £ d^ <-Cl < 2 ' 



In the second case, the right hand side of (|46|) becomes 

_ da df_ p2_df_ 

dx dx pi dx 



(48) 



Again, by choosing pi/vi = 1 + £, with e > 0, we have the following upper bound for the 
expression in ([35]): 

da df df 
dx dx dx 

Thus, it follows that fj,p t i (J) < — c 2 provided that the above quantity is uniformly negative 
definite. Since, by hypotheses, 

fin f) f 

(50) 



da n df u _ 
ox ox 



then 3c| 2 



"ff~'~§z~'~pr§z — — c 2 2- The proof of the Theorem is now complete. 



□ 
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From a biological viewpoint, the hardest hypothesis to satisfy in Theorem [6] might be that 
on the derivatives of / (x, y). However, it is possible to relax the hypothesis on df /dx if the rate 
of change of a (x) with respect to x, i.e. da/dx, is sufficiently larger than df/dx. In particular, 
the following result can be proved. 

Theorem 7. System without inputs u, evolving on a convex compact set, is contractive 

provided that: 

• da/dx > 0, Vx G C; 

. df/dy >0, Vx,y eC; 

• da/dx > maxc {2 \df /dx\} . 

Proof. The proof is similar to that of Theorem EJ In particular, we can repeat the same 
derivation to obtain again inequality (|46p . Thence, as no hypothesis is made on the sign of 
df/dx, choosing P2/P1 = 1 + £ we have 



da df 
h — + 

dx dx 



P2_df_ 

Pi dx 



da df 
1 — — + 

dx dx 



df 



dx 



+ e 



df 



dx 



(51) 



Thus, it follows that, if da/dx > 2\df/dx\, then 3 c 2 such that Hp,i(J) < — c 2 , implying 
contractivity. The above condition is satisfied by hypotheses, hence the theorem is proved. □ 

Remarks 

Theorems [6] and [7] show the possibility of designing with high flexibility the self-degradation 
and interaction functions for an input-output module. 

This flexibility can be further increased, for example in the following ways: 

• Results similar to that of the above Theorems can be derived (and also extended) if some 
self degradation rate for y is present in (jlOl) . i.e. 



x = u(t) -a(x) + f(x,y) 
V = -b(y) - f(x,y) 



(52) 



with g < 0. 



• Theorem [6] and Theorem [7] can also be extended to the case in which the X-module drives 
more than one downstream transcriptional modules. 

3 Materials and Methods 



All simulations are performed in MATLAB (Simulink), Version 7.4, with variable step ODE 
solver ODE23t. Simulink models are available upon request. 
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4 Conclusions 



We have presented a systematic methodology to derive conditions for transcriptional modules 
to be globally entrained to periodic inputs. By means of contraction theory, a useful tool from 
dynamical systems, we showed that it is possible to use non-Euclidean norms and their associ- 
ated matrix measures to characterize the behavior of several modules when subject to external 
periodic excitations. Specifically, starting with a simple bimolecular reaction, we considered 
the case of a general externally-driven transcriptional module and extended the analysis to 
some important generalizations including the case of multiple driven systems. In all cases con- 
ditions are derived by proving that the module of interest is contracting under some generic 
assumptions on its parameters. The importance of the results presented in the paper from a 
design viewpoint are also discussed by means of more abstract systems where generic nonlinear 
degradation and interaction terms are assumed. 
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A il-reachable sets 

We will make use of the following definition: 
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Definition 1. Let K > be any positive real number. A subset C C W 1 is fT-reachable if, for 
any two points xq and yo in C there is some continuously differentiable curve 7 : [0, 1] — ► C 
such that: 

1. 7 (0) =x , 

2. 7 (1) = y and 

3. \i (r) I <K\y Q -x Q \, Vr. 

For convex sets C, we may pick 7(7") = xo + r(yo — xq), so 7'(r) = yo — ^0 an d we can take 
K = 1. Thus, convex sets are 1-reachable, and it is easy to show that the converse holds as 
well. 

Notice that a set C is iT-reachable for some K if and only if the length of the geodesic 
(smooth) path (parametrized by arc length), connecting any two points x and y in C, is bounded 
by some multiple Kq of the Euclidean norm, \y — x\ 2 . Indeed, re-parametrizing to a path 7 
defined on [0, 1], we have: 

\l' (r)\ 2 < K \y-x\ 2 . 

Since in finite dimensional spaces all the norms are equivalent, then it is possible to obtain a 
suitable K for Definition [TJ 

Remark 1. The notion of K -reachable set is weaker that that of convex set. Nonetheless, in 
Theorem^ we will prove that trajectories of a smooth system, evolving on a K -reachable set, 
converge towards each other, even if C is not convex. This additional generality allows one to 
establish contracting behavior for systems evolving on phase spaces exhibiting "obstacles" , as 
are frequently encountered in path-planing problems, for example. A mathematical example of 
a set with obstacles follows. 

Example 1. Consider the two dimensional set, C , defined by the following contraints: 

x 2 + y 2 >l, x>0, y > . 

Clearly, C is a non-convex subset ofM?. We claim that C is K-reachable, for any positive real 
number K > ^ . Indeed, given any two points a and b in C , there are two possibilities: either 
the segment connecting a and b is in C , or it intersects the unit circle. In the first case, we can 
simply pick the segment as a curve (K = 1). In the second case, one can consider a straight 
segment that is modified by taking the shortest perimiter route around the circle; the length of 
the perimeter path is at most — times the length of the omitted segment. (In order to obtain a 
differentiable, instead of merely a piecewise- differentiable, path, an arbitrarily small increase in 
K is needed.) 

B Proof of Theorem [I] 

We now prove the main result on contracting systems, i.e. Theorem [U under the hypotheses 
that the set C, i.e. the set on which the system evolves, is .ff -reachable. 

Theorem 8. Suppose that C is a K-reachable subset ofM. n and that f(t,x) is infinitesimally 
contracting with contraction rate c 2 . Then, for every two solutions x(t) = (£>(i,0,£) and z(t) = 
ip(t, 0, () it holds that: 

\x(t) -z(t)\ < Ke- cH |f -C| Vi>0. (53) 
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Proof. Given any two points x (0) = £ and z (0) = £ in C, pick a smooth curve 7 : [0, 1] — ► C, 
such that 7 (0) = £ and 7 (1) = (. Let ?/> (i, r) = 0, 7 (r), that is, the solution of system ([T|) 
rooted in tp (0, r) = 7 (r), r e [0, 1]. Since and 7 are continuously differentiable, also ip (t,r) 
is continuously differentiable in both arguments. We define 

. , dip . . 
:= — (t,r). 



It follows that 



Now, 



so, we have: 



dw , . d f dip\ d { dip\ d „ , . , . . 



|-/»((,r),() = |^W(t,r),*)^((,r) 
or ox or 



-^(t,r) = J^(t,r),t)w(t,r), (54) 
where J(ip (t,r) ,t) = ||(^ (i, r) , i). Using Coppel's inequality (Vidyasagar 1993), yields 

\w(t,r)\ < |u»(0,r)|e^^ J(r))dT < K \£ - C\ e~ cH , (55) 
\/x E C, Vt € R + , and Vr € [0, 1]. Notice the Fundamental Theorem of Calculus, we can write 

ip(t, 1) - V(*,0) = y w(t,s)ds. 

Hence, we obtain 

-1 







|x(£) - < / |w(i,s)[ds. 
Jo 

Now, using ([55]) . the above inequality becomes: 

\x(t)-z(t)\ < J (|^(0,s)|e/oKJ(T))^ dg < K|e-C|e" c2 *. 
The Theorem is then proved. □ 



Proof of Theorem [TJ The proof follows trivially from Theorem El after having noticed that 
in the convex case, we may assume K = 1. □ 

C Proof of Theorem [2] 

In this Section we assume that the vector field / is T-periodic and prove Theorem [2j 

Before starting with the proof of Theorem [2] we make the following: 
Remark 2. Periodicity implies that the initial time is only relevant modulo T . More precisely: 

tp(kT + t,kT,£) = <p(t,0,£) Vk e N,t > 0, x G C. (56) 
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Indeed, let z(s) = (p(s, kT,£), s > kT, and consider the function x(t) = z{kT + t) = <p(kT + 
t,kT,£), fort > 0. So, 

±{t) = z(kT + t) = f(kT + t, z(kT + t)) = f(kT + t, x(t)) = f(t, x{t)) , 

where the last equality follows by T -periodicity of f . Since x(0) = z(kT) = <p(kT, kT,£) = £, it 
follows by uniqueness of solutions that x(t) = ip(t,0,£) = <p(kT + t, AT, £), which is [56\). As a 
corollary, we also have that 

<p(kT + 1, 0, f ) = ip(kT + 1, kT, <p(kT, 0, £)) = <p(t, 0, <p(kT, 0,0) V k e N, t > 0, x £ C (57) 

where the first equality follows from the semigroup property of solutions (see e.g. (Sontag 1998)), 
and the second one from $56)) applied to ip(kT,0,£) instead of £. 

Define now 

p(0 = <p(t,o,€), 

where £ = x (0) G C. The following Lemma will be useful in what follows. 
Lemma 1. P k (£) = <p(kT,0,€) for allkeN and £ G C. 

Proof. We will prove the Lemma by recursion. In particular, the statement is true by definition 
when k = 1. Inductively, assuming it true for k, we have: 

P fc+1 (£) = P(P k (0) = <P(T, 0, P k (0) = <P(T, 0, ip(kT, 0, 0) = v(kT + T, 0, , 
as wanted. □ 
Theorem 9. Suppose that: 

• C is a closed K -reachable subset o/R"; 

m f is infinitesimally contracting with contraction rate c 2 ; 

• f is T -periodic; 
. Ke~ c2T < 1. 

Then, there is an unique periodic solution a(t) : [0, oo) — > C of fljj having period T. Further- 
more, every solution x{t), such that x (0) = £ € C , converges to a (t), i.e. \x(t) — a(t)\ — > as 
t — ► oo. 

Proof. Observe that P is a contraction with factor Ke~ c2T < 1: |P(£) — P(C)I ^ Ke~ c2T \£ — £| 
for all £, ( G C, as a consequence of Theorem [SJ The set C is a closed subset of W 2 and 
hence complete as a metric space with respect to the distance induced by the norm being 
considered. Thus, by the contraction mapping theorem, there is a (unique) fixed point £ of 
P. Let a(t) := (p(t,Q,£). Since a(T) = P(q) = £ = a(0), a(t) is a periodic orbit of period T. 
Moreover, again by Theorem [HJ we have that \x(t) — a(t)\ < Ke~ c * £ — £ — > 0. Uniqueness 
is clear, since two different periodic orbits would be disjoint compact subsets, and hence at 
positive distance from each other, contradicting convergence. This completes the proof. □ 
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Proof of Theorem [2} It will suffice to note that the assumption Ke c T < 1 in Theorem [9] is 
automatically satisfied when the set C is convex (i.e. K = 1) and the system is infinitesimally 



Notice that, even in the non-convex case, the assumption Ke~ c < 1 can be ignored, if 
we are willing to assert only the existence (and global convergence to) a unique periodic orbit, 
with some period kT for some integer k > 1. Indeed, the vector field is also fcT-periodic for 
any integer k. Picking k large enough so that Ke~ c kT < 1, we have the conclusion that such 
an orbit exists, applying Theorem [9l 

D Cascades 

In order to show that cascades of contracting systems remain contracting, it is enough to show 
this, inductively, for a cascade of two systems. 

Consider a system of the following form: 



where x(t) G C\ C f" 1 and y(t) G C 2 C W 12 for all t (C% and C 2 are two K-reachable sets). We 
write the Jacobian of / with respect to x as A(t, x) = x), the Jacobian of g with respect to 
x as B(t,x,y) = ^(t,x,y), and the Jacobian of g with respect to y as C(t,x,y) = ||(i,x,y), 

We assume the following: 

1. The system x = f(t,x) is infinitesimally contracting with respect to some norm (generally 
indicated as |»|*), with some contraction rate cf, that is, fi*(A(t, x)) < — cf for all x € C\ 
and all t > 0, where /i* is the matrix measure associated to \*\*- 

2. The system y = f(t,x,y) is infinitesimally contracting with respect to some norm (which 
is, in general different from l*^, and is denoted by !•!„), with contraction rate c|, when x 
is viewed a a parameter in the second system, that is, fi^(C(t, x, y)) < —c\ for all x E C\, 
y G C2 and all t > 0, where /i** is the matrix measure associated to 1*1^. 

3. The mixed Jacobian B(t, x, y) is bounded: \\B(t, x, y) \\ < k 2 , for all x E C\, y G C2 and all 
i > 0, for some real number k, where "||»||" is the operator norm induced by 1*1^ and |«|^ 
on linear operators W llXri2Xl — > R niX ™ 2 . (All norms in Euclidean space being equivalent, 
this can be verified in any norm.) 

We claim that, under these assumptions, the complete system is infinitesimally contracting. 
More precisely, pick any two positive numbers p\ and p 2 such that 



contracting. 



□ 



x = f(t,x) 

y = g(t,x,y) 




and let 
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We will show that //(J) < — c 2 , where J is the full Jacobian: 

J 



A 

B C 



with respect to the matrix measure \i induced by the following norm in M. niXU2 : 

\(xi,x 2 )\ =Pi +P2 \x 2 \^ ■ 

Since 

(j + hJ ) x -\ ( J + hA ) x * 
1 + ' I hB Xl + (/ + hC)x 2 

for all h and x, we have that, for all h and x: 

\{I + hJ)x\ = p 1 \(I + hA)x 1 \+p 2 \hBx 1 + (I + hC)x 2 \ 

< pi \I + |xi| +p 2 \hB\ \xi \ +p 2 \I + /iC| |x 2 | , 

where from now on we drop subscripts for norms. Pick now any h > and a unit vector x 
(which depends on h) such that \\I + hJ\\ = \(I + hJ)x\. Such a vector x exists by the definition 
of induced matrix norm, and we note that 1 = \x\ = p\ \x\\^ +p 2 \x 2 \^, by the definition of the 
norm in the product space. Therefore: 



\{\\I + hJ\\ 



= ±{\{I + hJ)x\-\x\) 

< ^ (Pi \I + hA\ \xi \ +p 2 \hB\ \x\ \ +p 2 \I + hC\ \x 2 \ —p\ \x\\ —p 2 \x 2 \) 
\ (\I + hA\-l + ^h \B\) Pl \ Xl \ + \ (|J + hC\ - l)p2 \x 2 \ 



h 



Pi 



h 



< max<i y(\I + hA\ + — k 2 , \ (\I + hC\ 
1 h pi h 



1) 



where the last inequality is a consequence of the fact that Aiai + \ 2 a2 < max{ai,a2} for any 
nonnegative numbers with Ai + A2 = 1 (convex combination of the a^s). Now taking limits as 
h \ 0, we conclude that 



fi( J) < max < — c? + — k 2 , — Cn \ 
I Pi J 



as desired. 



E A counterexample to entrainment 

In (Sontag 2009) there is given an example of a system with the following property: when the 
external signal u(t) is constant, all solutions converge to a steady state; however, when u(t) = 
s'mt, solutions become chaotic. (Obviously, this system is not contracting.) The equations are 
as follows: 

x = —x — u 

p = — p + a(x + u) 

i = io(V-£) 

ip = 28K - V> - p£( 

C = P^-(8/3)C 
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where a(y) = y 2 j(K + y 2 ) and K = 0.0001. Figure [8] shows typical solutions of this system 
with a periodic and constant input respectively. The function "rand" was used in MATLAB to 
produce random values in the range [—10, 10]. 




10 20 30 40 50 60 70 80 90 100 10 20 30 40 50 60 70 80 90 100 



Figure 8: Simulation of counter-example, done with the following randomly-chosen input and 
initial conditions: u{t) = 1.89, x(0) = 2.95 p(0) = -0.98, f (0) = 0.94, ^(0) = -4.07, ((0) = 
4.89. Green: inputs are u(t) = sint (left panel) and u(t) = 5.13 (randomly picked, right panel). 
Blue: Note chaotic-like behavior in response to periodic input, but steady state in response 
to constant input. 
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